extern "C" {
#include "cblas.h"
}
#include "ParallelMatrixMultiply.h"
#include "Matrix.h"
#include <omp.h>
#include <iostream>
#include <math.h>

using namespace std;

MatrixMultiply::MatrixMultiply(Matrix const *matA, Matrix const *matB, Matrix *matC) : m_matA(*matA), m_matB(*matB), m_matC(*matC), n(matA->GetRowCount()), m_lda(matA->GetLD()), m_ldb(matB->GetLD()), m_ldc(matC->GetLD())
{
	m_pMatA = matA->GetMatrix();
	m_pMatB = matB->GetMatrix();
	m_pMatC = matC->GetMatrix();
}

void MatrixMultiply::method1(const int threads)
{
	int n = m_matA.GetRowCount();
	int n1=n/2;
	int n2=n/2;
	int m1=n/2;
	int m2=n/2;
	int k1=n/2;
	int k2=n/2;	
	if(n%2)
	{
		n1=int(n/2)+1;
		m1=int(n/2)+1;
		k1=int(n/2)+1;
	}
	double alpha =1.0;
	double beta0 =0.0;
	double beta1 =1.0;
	//std::cout << "n: " << n << "; Rechnung: " << int(n/2)+1 << std::endl;
	//std::cout << "n1: " << n1 << "; n2: " << n2 << std::endl;
	//std::cout << "m1: " << m1 << "; m2: " << m2 << std::endl;
	//std::cout << "k1: " << k1 << "; k2: " << k2 << std::endl;
	//std::cout << std::endl;
	#pragma omp parallel num_threads(threads) private(n)
	{
		int indexA1;
		int indexA2;
		int indexB1;
		int indexB2;
		int indexC;
		int m;
		int k;
		#pragma omp sections
		{
			#pragma omp section
			{
				//left partial matrix
				//Berechnung C11
				indexA1 = 0;
				m = m1;
				indexB1 = 0;
				k = k1;
				indexC = 0;
				n = n1;
				//std::cout << "C11 Matrix A11 Soll(1) Ist Ptr: " << m_pMatA[indexA1] << std::endl;
				//std::cout << "C11 Matrix B11 Soll(1) Ist Ptr: " << m_pMatB[indexB1] << std::endl;
				//std::cout << "C11 Matrix C Ptr: " << m_pMatC[indexC] << std::endl;
			
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, &m_pMatA[indexA1], m_lda, &m_pMatB[indexB1], m_ldb, beta0, &m_pMatC[indexC], m_ldc);
				//right  partial matrix
				indexA2 = m_lda*k1;
				m = m1;
				indexB2 = k1;
				k = k2;
				indexC = 0;
				n = n1;
				//std::cout << "C11 Matrix A12 Soll(16) Ist Ptr: " << m_pMatA[indexA2] << std::endl;
				//std::cout << "C11 Matrix B21 Soll(4) Ist Ptr: " << m_pMatB[indexB2] << std::endl;
				//std::cout << "C11 Matrix C Ptr: " << m_pMatC[indexC] << std::endl;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k, alpha, &m_pMatA[indexA2], m_lda, &m_pMatB[indexB2], m_ldb, beta1, &m_pMatC[indexC], m_ldc);
			}
			#pragma omp section
			{
				//Berechnung C12
				indexA1 = 0;
				m = m1;
				indexB1 = n1*m_ldb;
				k = k1;
				indexC = n1*m_ldc;
				n = n2;
				//std::cout << "C12 Matrix A11 Soll(1) Ist Ptr: " << m_pMatA[indexA1] << std::endl;
				//std::cout << "C12 Matrix B12 Soll(16) Ist Ptr: " << m_pMatB[indexB1] << std::endl;
				//std::cout << "C12 Matrix C Ptr: " << m_pMatC[indexC] << std::endl;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k , alpha, &m_pMatA[indexA1], m_lda, &m_pMatB[indexB1], m_ldb, beta0, &m_pMatC[indexC], m_ldc);
				//right  partial matrix
				indexA2 = m_lda*k1;
				m = m1;
				indexB2 = (n1*m_ldb)+k1;
				k = k2;
				indexC = n1*m_ldc;
				n = n2;
				//std::cout << "C12 Matrix A12 Soll(16) Ist Ptr: " << m_pMatA[indexA2] << std::endl;
				//std::cout << "C12 Matrix B22 Soll(19) Ist Ptr: " << m_pMatB[indexB2] << std::endl;
				//std::cout << "C12 Matrix C Ptr: " << m_pMatC[indexC] << std::endl;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k, alpha, &m_pMatA[indexA2], m_lda, &m_pMatB[indexB2], m_ldb, beta1, &m_pMatC[indexC], m_ldc);
			}
			#pragma omp section
			{
				//Berechnung C21
				indexA1 = m1;
				m = m2;
				indexB1 = 0;
				k = k1;
				indexC = m1;
				n = n1;
				//std::cout << "C21 Matrix A21 Soll(4) Ist Ptr: " << m_pMatA[indexA1] << std::endl;
				//std::cout << "C21 Matrix B11 Soll(1) Ist Ptr: " << m_pMatB[indexB1] << std::endl;
				//std::cout << "C21 Matrix C Ptr: " << m_pMatC[indexC] << std::endl;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k , alpha, &m_pMatA[indexA1], m_lda, &m_pMatB[indexB1], m_ldb, beta0, &m_pMatC[indexC], m_ldc);
				//right  partial matrix
				indexA2 = m_lda*k1+m1;
				m = m2;
				indexB2 = k1;
				k = k2;
				indexC = m1;
				n = n1;
				//std::cout << "C21 Matrix A22 Soll(19) Ist Ptr: " << m_pMatA[indexA2] << std::endl;
				//std::cout << "C21 Matrix B21 Soll(4) Ist Ptr: " << m_pMatB[indexB2] << std::endl;
				//std::cout << "C21 Matrix C Ptr: " << m_pMatC[indexC] << std::endl;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k, alpha, &m_pMatA[indexA2], m_lda, &m_pMatB[indexB2], m_ldb, beta1, &m_pMatC[indexC], m_ldc);
			}
			#pragma omp section
			{
				//Berechnung C22
				indexA1 = m1;
				m = m2;
				indexB1 = n1*m_ldb;
				k = k1;
				indexC = n1*m_ldc+m1;
				n = n2;
				//std::cout << "C22 Matrix A21 Soll(4) Ist Ptr: " << m_pMatA[indexA1] << std::endl;
				//std::cout << "C22 Matrix B12 Soll(16) Ist Ptr: " << m_pMatB[indexB1] << std::endl;
				//std::cout << "C22 Matrix C Ptr: " << m_pMatC[indexC] << std::endl;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k , alpha, &m_pMatA[indexA1], m_lda, &m_pMatB[indexB1], m_ldb, beta0, &m_pMatC[indexC], m_ldc);
				//right  partial matrix
				indexA2 = (m_lda*k1)+m1;
				m = m2;
				indexB2 = (n1*m_ldb)+k1;
				k = k2;
				indexC = n1*m_ldc+m1;
				n = n2;
				//std::cout << "C22 Matrix A22 Soll(19) Ist Ptr: " << m_pMatA[indexA2] << std::endl;
				//std::cout << "C22 Matrix B22 Soll(19) Ist Ptr: " << m_pMatB[indexB2] << std::endl;
				//std::cout << "C22 Matrix C Ptr: " << m_pMatC[indexC] << std::endl;
				cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n ,k, alpha, &m_pMatA[indexA2], m_lda, &m_pMatB[indexB2], m_ldb, beta1, &m_pMatC[indexC], m_ldc);
			}
		}
	}
}

int workload(int problem_size, int num_threads, int thread_id)
{
	int basic = floor(problem_size/num_threads);
	int rest = (problem_size%num_threads);
	basic += thread_id<rest;
	return basic;
}

int workload_startaddress(int problem_size, int num_threads, int thread_id)
{
	int basic = floor(problem_size/num_threads);
	int rest = max(0, (problem_size%num_threads));
	return min(rest, thread_id)*(basic+1)+max(0,thread_id-rest)*basic;
}

void MatrixMultiply::method2(const int threads)
{
	int m = m_matB.GetColumnCount(); // must be equal to m_matA/C.GetRowCount();
	int n = m_matA.GetRowCount(); // must be equal to m_matB/C.GetColumnCount();
	int k = m_matA.GetColumnCount(); // must be equal to m_matB/C.GetRowCount(); and m_matA.GetColumnCount();

	int indexA = 0;
	int num_threads = threads;
	int problem_size = m_matA.GetRowCount();
	int *problems = new int[num_threads];
	int *starts = new int[num_threads];

	double alpha =1.0;
	//double beta0 =0.0;
	double beta1 =0.0;

	for(int i=0 ; i < num_threads ; ++i)
	{
		problems[i] = workload(problem_size, num_threads, i);
		starts[i] = m_ldb*workload_startaddress(problem_size, num_threads, i);
	}
//	cout << "probs: ";
//	for(int i=0 ; i < num_threads ; ++i)
//	{
//		cout << problems[i] << " ";
//	}
//	cout << endl <<"starts: ";
//	for(int i=0 ; i < num_threads ; ++i)
//	{
//		cout << starts[i] << " ";
//	}
	#pragma omp parallel for num_threads(num_threads) shared(problems, starts)
//	#pragma omp parallel for shared(problems, starts)
	for ( int i = 0 ; i < num_threads ; ++i)
	{
		int indexBC= starts[i];
		cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, problems[i], k, alpha, &m_pMatA[indexA], m_lda, &m_pMatB[indexBC], m_ldb, beta1, &m_pMatC[indexBC], m_ldc);
	}

	//double alpha = 1.0;
	//double beta = 0.0;
	//const int n = m_matA.GetRowCount();
	//cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, alpha, &m_pMatA[0], m_lda, &m_pMatB[0], m_ldb, beta, &m_pMatC[0], m_ldc);
	delete [] problems;
	delete [] starts;
}

void MatrixMultiply::method3(const int threads)
{
	int k = m_matA.GetColumnCount();

	int num_threads = threads;
	int problem_sizeCol = m_matB.GetColumnCount();	//n
	int problem_sizeRow = m_matA.GetRowCount();	//m
	int *problemsCol = new int[num_threads];
	int *startsCol = new int[num_threads];
	int *problemsRow = new int[num_threads];
	int *startsRow = new int[num_threads];

	double alpha =1.0;
	//double beta0 =0.0;
	double beta1 =0.0;

	for(int i=0 ; i<num_threads ; ++i)
	{
		problemsCol[i] = workload(problem_sizeCol, num_threads, i);
		startsCol[i] = workload_startaddress(problem_sizeCol, num_threads, i);
		problemsRow[i] = workload(problem_sizeRow, num_threads, i);
		startsRow[i] = m_ldb*workload_startaddress(problem_sizeRow, num_threads, i);
	}
//	cout << endl << "probs: ";
//	for(int i=0 ; i < num_threads ; ++i)
//	{
//		cout << problemsCol[i] << " ";
//	}
//	cout << "starts: ";
//	for(int i=0 ; i < num_threads ; ++i)
//	{
//		cout << startsCol[i] << " ";
//	}
	for(int a=0 ; a<num_threads ; ++a)
	{
		#pragma omp parallel for num_threads(num_threads) shared(problemsCol, startsCol, problemsRow, startsRow)
		for ( int i=0 ; i<num_threads ; ++i)
		{
			int j = (a+i)%num_threads;
			int indexC = startsRow[j]+startsCol[i];
			int indexA = startsCol[i];
			int indexB = startsRow[j];
			int m = problemsCol[i];
			int n = problemsRow[j];
			//int k = k
			cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, &m_pMatA[indexA], m_lda, &m_pMatB[indexB], m_ldb, beta1, &m_pMatC[indexC], m_ldc);
		}
	}
	delete [] problemsCol;
	delete [] startsCol;
	delete [] problemsRow;
	delete [] startsRow;
}

void MatrixMultiply::cblas_myDGEMM(const int threads)
{
	double alpha = 1.0;
	double beta = 0.0;
	const int n = m_matA.GetRowCount();
	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, alpha, &m_pMatA[0], m_lda, &m_pMatB[0], m_ldb, beta, &m_pMatC[0], m_ldc);
}

MatrixMultiply::~MatrixMultiply(void)
{
}
